### Hollibaugh and Rothenberg
### Agency Control through the Appointed Hierarchy: Presidential Politicization of Unilateral Appointees
### Utility Functions


### utility function --- create texreg table for POLR objects

extract.polr <- function(model, include.thresholds = FALSE, include.aic = TRUE,
                         include.bic = TRUE, include.loglik = TRUE,
                         include.deviance = FALSE,
                         include.nobs = TRUE, include.lrtest = TRUE, ...) {
  s <- summary(model, ...)
  lrtest.chi <- round(lmtest::lrtest(model)$C[2], 3)
  lrtest.pval <- round(lmtest::lrtest(model)$P[2], 3)

  tab <- s$coefficients
  zeta.names <- names(s$zeta)
  beta <- tab[!rownames(tab) %in% zeta.names, ]
  thresh <- tab[rownames(tab) %in% zeta.names, ]

  if (sum(!rownames(tab) %in% zeta.names) == 1) {
    beta <- t(beta)
    rownames(beta) <- rownames(tab)[!rownames(tab) %in% zeta.names]
  }
  if (sum(rownames(tab) %in% zeta.names) == 1) {
    thresh <- t(thresh)
    rownames(thresh) <- rownames(tab)[rownames(tab) %in% zeta.names]
  }

  threshold.names <- rownames(thresh)
  threshold.coef <- thresh[, 1]
  threshold.se <- thresh[, 2]
  threshold.zval <- thresh[, 1] / thresh[, 2]
  threshold.pval <- 2 * pnorm(abs(threshold.zval), lower.tail = FALSE)

  beta.names <- rownames(beta)
  beta.coef <- beta[, 1]
  beta.se <- beta[, 2]
  beta.zval <- beta[, 1] / beta[, 2]
  beta.pval <- 2 * pnorm(abs(beta.zval), lower.tail = FALSE)

  if (include.thresholds == TRUE) {
    names <- c(beta.names, threshold.names)
    coef <- c(beta.coef, threshold.coef)
    se <- c(beta.se, threshold.se)
    pval <- c(beta.pval, threshold.pval)
  } else {
    names <- beta.names
    coef <- beta.coef
    se <- beta.se
    pval <- beta.pval
  }

  n <- nobs(model)
  lik <- logLik(model)[1]
  aic <- AIC(model)
  bic <- BIC(model)
  dev <- deviance(model)
  lrt <- lrtest.chi
  lrp <- lrtest.pval

  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  if (include.aic == TRUE) {
    gof <- c(gof, aic)
    gof.names <- c(gof.names, "AIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.bic == TRUE) {
    gof <- c(gof, bic)
    gof.names <- c(gof.names, "BIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.loglik == TRUE) {
    gof <- c(gof, lik)
    gof.names <- c(gof.names, "Log Likelihood")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.deviance == TRUE) {
    gof <- c(gof, dev)
    gof.names <- c(gof.names, "Deviance")
    gof.decimal <- c(gof.decimal, TRUE)
  }

  if (include.lrtest == TRUE) {
    gof <- c(gof, lrt, lrp)
    gof.names <- c(gof.names, "Likelihood Ratio Test", "%Likelihood Ratio p-value")
    gof.decimal <- c(gof.decimal, TRUE, TRUE)
  }

  if (include.nobs == TRUE) {
    gof <- c(gof, n)
    gof.names <- c(gof.names, "Number of Observations")
    gof.decimal <- c(gof.decimal, FALSE)
  }

  tr <- createTexreg(
    coef.names = names,
    coef = coef,
    se = se,
    pvalues = pval,
    gof.names = gof.names,
    gof = gof,
    gof.decimal = gof.decimal
  )
  return(tr)
}

setMethod("extract", signature = className("polr", "MASS"),
          definition = extract.polr)


# get Quasi-AIC for quasibinomial models
x.quasibinomial <- function(...) {
  res <- quasibinomial(...)
  res$aic <- binomial(...)$aic
  res
}


### utility function --- quasi-AIC for quasiGLM objects
manual_QAIC <- function(x){QAIC(update(x, family = x.quasibinomial),
                                chat = deviance(update(x, family = x.quasibinomial))/
                                  df.residual(update(x, family = x.quasibinomial)))}


### utility function --- create texreg table for GLM objects
# extension for glm objects
extract.glm <- function(model, include.aic = FALSE, include.bic = FALSE,
                        include.loglik = FALSE, include.deviance = FALSE,
                        include.qaic = TRUE,
                        include.wald = TRUE,
                        include.nobs = TRUE, ...) {
  s <- summary(model, ...)

  coefficient.names <- rownames(s$coef)
  coefficients <- s$coef[, 1]
  standard.errors <- s$coef[, 2]
  significance <- s$coef[, 4]

  aic <- AIC(model)
  bic <- BIC(model)
  qaic <- manual_QAIC(model)
  waldstat <- lmtest::waldtest(model, test = "Chisq")$C[2]
  waldpval <- lmtest::waldtest(model, test = "Chisq")$P[2]
  lik <- logLik(model)[1]
  dev <- deviance(model)
  n <- nobs(model)

  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  if (include.aic == TRUE) {
    gof <- c(gof, aic)
    gof.names <- c(gof.names, "AIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.qaic == TRUE) {
    gof <- c(gof, qaic)
    gof.names <- c(gof.names, "Quasi-AIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.bic == TRUE) {
    gof <- c(gof, bic)
    gof.names <- c(gof.names, "BIC")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.loglik == TRUE) {
    gof <- c(gof, lik)
    gof.names <- c(gof.names, "Log Likelihood")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.deviance == TRUE) {
    gof <- c(gof, dev)
    gof.names <- c(gof.names, "Deviance")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.wald == TRUE) {
    gof <- c(gof, waldstat, waldpval)
    gof.names <- c(gof.names, "Wald Test", "%Wald Test p-value")
    gof.decimal <- c(gof.decimal, TRUE, TRUE)
  }
  if (include.nobs == TRUE) {
    gof <- c(gof, n)
    gof.names <- c(gof.names, "Number of Observations")
    gof.decimal <- c(gof.decimal, FALSE)
  }

  tr <- createTexreg(
    coef.names = coefficient.names,
    coef = coefficients,
    se = standard.errors,
    pvalues = significance,
    gof.names = gof.names,
    gof = gof,
    gof.decimal = gof.decimal
  )
  return(tr)
}

setMethod("extract", signature = className("glm", "stats"),
          definition = extract.glm)


